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ABSTRACT 


This thesis examines two implementations of the para- 
bolic equation approximation to the acoustic wave equation 
aimed at removing three errors inherent to the wide-angle 
parabolic equation (WAPE) model. First, the selection of 
the range-step size used by the split-step Fourier algorithm 
affects the convergence of the solution. Second, in certain 
ocean environments WAPE incorrectly computes the down-range 
transmission loss. Finally, WAPE does not reproduce the 
standard normal mode basis set as defined by normal mode the- 
ory. A double-precision implementation of the WAPE (DP- 
WAPE) is developed to evaluate the dependence of solution 
convergence on the numerical precision of the model. 
Finally, an implementation that is insensitive to the choice 
of the reference sound speed (COIPE) is evaluated for its 
ability to reduce or remove the latter two of these three 
errors. The stability of the WAPE solution was found to be 
unaffected by the DP-WAPE implementation. The range-step 
dependence is inherent to the split-step algorithm. The 
COIPE corrects the transmission loss anomaly and satisfacto- 


rily reproduces the standard normal mode basis set. 
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I. INTRODUCTION 


The fall of the former Soviet Union and the resulting 
shift in the international political power base has resulted 
in a dramatic shift in both the content and the application 
of United States Naval forces during the past decade. The 
world’s largest blue water navy finds itself increasingly 
involved in green water operations. This siritsto littoral 
waters has brought to the forefront the ever present danger 
of easily obtainable and extremely capable diesel submarines 
and shallow water mines. Countries with otherwise limited 
military resources can acquire the ability to thwart or at 
the least hinder the application of national objectives. As 
a result of these changes, numerous efforts toward the 
development of more accurate acoustic propagation prediction 
models have arisen. These include areas such as matched- 
field processing, transient localization, and underwater 
acoustic communications. Essential to these efforts is the 
development of acoustic prediction models which can produce 
both valid and accurate results in highly variable shallow 


water ocean environments. 


A. OCEAN ACOUSTIC MODELING TECHNIQUES 


Current acoustic modeling efforts generally fall into 
one of four broad categories. These are ray models, wave- 
number integration techniques, normal mode models, and para- 
bolic equation models. Each of these are based on a variety 
of approximations to the linear acoustic wave equation. 

Many of the early models are based on the geometrical limit 
and define ray trajectories of sound propagation. The 
resulting ray models can quickly predict pulse propagation 
travel times and provide rough approximations of the sound 
field. The limitation of this approach is its reliance ona 
high-frequency limit which neglects finite-frequency effects 
such as diffraction. In addition, it has been shown (Smith, 
et. al., 1992) that the ray equations form a set of coupled, 
nonlinear equations which suffer from chaotic solutions in 
the presence of range-dependence. Given the current impor- 
tance of low-frequency propagation in range-dependent envi- 
ronments, the usefulness of ray models in shallow water 
media is limited. 

Wavenumber integration methods produce highly accurate, 


full-wave solutions to the wave equation. However, they are 


computationally intensive and are, by design, limited to 
range-independent environments. They can be quite useful in 
determining near-field effects but are not of practical 
utility in predictions of far-field acoustic propagation. 

Normal mode models are based on a separation of vari- 
ables in range and depth of the acoustic wave equation. 
This approach maintains the full-wave characteristics of the 
acoustic field, and provides highly accurate results to both 
simple and complex range dependent ocean environments. The 
solutions produced are based on the single-frequency Helm- 
holtz wave equation. Due to the computational complexity 
involved with high frequencies and broadband pulses, as well 
as range-dependent media, such models are best suited for 
low-frequency, shallow water environments which support only 
a small number of propagating modes. The results remain very 
accurate, but computational time can become counter produc- 
tive. 

The final acoustic modeling technique, and undoubtedly 
the most popular in range-dependent environments, is based 
on the parabolic approximation to the Helmholtz wave equa- 
tion. Numerous marching algorithms have been developed that 


solve the acoustic field as a boundary value problem. While 


) 


it 1S a one-way wave equation treatment, it is a full-wave 
model and therefore includes all of the finite-frequency 
effects of diffraction. Two-way adaptations of parabolic- 
equation models have been developed, but do not result in 
Significant changes for the majority of shallow water envi- 
ronments of interest. Solutions are, however, still single 
frequency, requiring multiple solution generation for pulse 
propagation responses. It is this method that will be the 


primary focus of this thesis. 


B. HISTORY OF THE PARABOLIC EQUATION APPROXIMATION 


The parabolic equation (PE) approximation method was 
first applied to the problem of underwater acoustics by Tap- 
pert (1977). Its historical roots are much deeper, however. 
Breakthroughs in physics and mathematics in the mid-1920’s 
provided the basis for the parabolic wave equation used in 
acoustics today. The standard parabolic equation (SPE) has 
the same form as the Schrodinger’s equation in quantum 
mechanics. The earliest documented use of the SPE as an 
approximation to the theory of wave propagation dates to the 
mid-1940’s and the work of Leontovich and Fock (1946) who 


originally coined the term “parabolic equation method”. 
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They applied the method to predicting the diffraction on 
long range tropospheric wave propagation caused by the 
spherical shape of the earth. Fock (1965) expanded this 
approach to include high-frequency scattering as well as 
microwave propagation in waveguides. 

With the advent of lasers and their coherent radiation 
sources in the 1960’s, the PE method was further extended to 
laser beam propagation where it is generally referred to as 
the “quasi-optical” equation. The quasi-optical equation is 
most widely used in nonlinear optics where the index of 
refraction depends on intensity. It is often called the non- 
linear Schrodinger’s equation. The PE method has also been 
applied, in more recent years, to the field of fiber optics 
and plasma physics. 

Beam propagation in random media also lends itself well 
to the PE method. Regardless of the source of the beam -- 
radio, acoustic, optical -- the problem is analogous to the 
quantum mechanics problem of motion in a random potential. 
The radar application of this problem in a randomly fluctu- 
ating ionosphere using the split-step Fourier algorithm, 
discussed later in the next chapter, is covered extensively 


by Hardin and Tappert (1974). 
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The most extensive use of the PE approximation to the 
elliptical wave equation has been in the field of low-fre- 
quency underwater acoustic propagation. More than 120 arti- 
cles and technical reports regarding the PE method and its 
applications to underwater acoustics were published in the 
last 15 years. Numerous marching algorithms and computer 
applications have been constructed with varying degrees of 
success. The earliest of these were based on the split-step 
Fourier algorithm (Hardin and Tappert, 1974). In general, 
these models accept input defining the sound speed, volume 
loss profiles, and depth contours, and produce output of the 
acoustic field and associated transmission loss (TL) curves. 
These results are generally in excellent agreement with 
experimental measurements for deep ocean environments. The 
theory behind this method and subsequent attempts to improve 
it are covered in detail in the next chapter. Since its ini- 
tial use, the PE method has been extended to include higher 
acoustic frequencies, random internal-wave fluctuations in 
the index of refraction, and shallow-water environments with 


mixed levels of success. 


Cc. ADVANTAGES AND LIMITATIONS OF THE PE METHOD 


Long-range, low-frequency sound propagation is gener- 
ally dominated by rays having small grazing angles since 
rays propagating at steep angles are greatly attenuated due 
to penetration and absorption in the seabed. Tappert intro- 
duced the PE method, which decomposes an elliptical wave 
equation into two equations through the separation of outgo- 
ing and incoming propagating fields. Tappert’s implementa- 
tion of this method was the first applied to underwater 
acoustics and is referred to as the standard parabolic equa- 
tion (SPE). The SPE implementation and many of its more 
recent incarnations share a number of inherent advantages 
and limitations. 

The advantages of the PE approximation can be divided 
into three categories (McDaniel, 1975). First, in the long- 
range, low-frequency environment outlined above, the govern- 
ing equation can be easily expressed by a parabolic equa- 
tion, which is easier to solve than either the elliptical or 
the hyperbolic types. The PE method, therefore, provides a 
relatively quick solution to long-range, low-frequency, 


range-dependent problems. Second, in solving the Helmholtz 


equation, the problem is posed as a pure boundary value prob- 
lem. In practice, the determination of the vertical bound- 
ary condition is often difficult and imprecise. The PE, 
however, poses the problem as one of an initial-boundary 
problem, therefore avoiding the difficulty of the vertical 
far-end wall boundary condition. Finally, the PE solution 
is generally cheaper to obtain in both memory and run-time 
than solving the elliptical equation. 

The limitations of the PE method can be placed in one of 
two broad types. The first of these results from the mathe- 
matical formulation of the PE itself which does not allow the 
inclusion of certain physical phenomenon. The PE is only 
valid in the far-field and therefore cannot be used for a 
complete analysis of near-field effects. The index of 
refraction is assumed to be slowly varying in range. Large 
and or abrupt changes in the index of refraction can limit 
the accuracy and validity of results. The PE uses the one- 
way wave equation and wave propagation from the reverse 
direction is ignored. Backscattering effects are, there- 
fore, neglected. Treatment of the so-called “Square root 
operator” defines the size of the propagation angle which in 


turn effects the validity of the PE. 
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The second type of limitation is defined by the method 
of solution to the PE. Many of the limitations of this type 
are specific to the solution algorithm being implemented. 
The fast Fourier transform works well only for equations 
with constant coefficients, and therefore must be used with 
caution. If the boundary is rigid, the algorithm cannot 
treat it exactly. Finite difference techniques require reg- 
ular grid partitions. Using the horizontal-interface PE to 
handle irregular interfaces, unless the sloping angle is 
small or density change is small, may not produce satisfac- 
tory results. In theory, there are no limitations on fre- 
quency. In practice, however, using finite alreerence 
techniques for high-frequency problems results in intolera- 
ble and expensive computational times. 

The two types of limitations result in three basic 
errors. The first error is due to the initial PE approxima- 
tion. Findings by Smith and Smith (1997) and others show 
that in the SPE a normal mode will be propagated with the 
correct amplitude and ee shape but with an error in phase 
velocity. Some higher order PE approximations improve the 
accuracy of the phase velocity, but may not properly propa- 


gate the mode amplitudes. The error in phase velocity is 


=) 


addressed in a number of PE approximation methods in the way 
in which the square root operator is expressed. The second 
error is introduced by the selection of the range step size. 
The error is typically found in the second or third order of 
the range step increment (Jensen, et. al., 1994). If the 


Geurer isg@efined as HE, then for the@ease*of severalmmedes 


propagating, the error produced in a single range step size 


would be E= SE,. The third type of error is due to trun- 


n 


cating the field. The initial approximation of the PE poses 
basic limitations as outlined above. The errors introduced 
by truncating the field are rooted mainly in computational 
efficiency in which the numerical results are obtained. A 
detailed analysis of the effects of range step choice were 


investigated by Smith (2000). 
D. THESIS SUMMARY 


The objective of this thesis is to explore the accuracy 
and validity of various approaches to the PE approximation 
method as it applies to underwater acoustics, in particular, 


how each of these methods or implementations affects the 
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errors outlined above. The resulting accuracy and validity 
of the model output will be analyzed. 

Chapter II, Theoretical Background, reviews the theory 
relevant to the material in the subsequent chapters. Over- 
views are provided of the parabolic equation approximation 
used for acoustic propagation modeling, the standard para- 
bolic equation (SPE) method, the wide angle parabolic equa- 


Elon method (WAPE);, and the relerence Soundware dm c.) 


insensitive parabolic equation method (COIPE) implementa- 
tions. Additionally, a brief overview of the method of nor- 
mal modes is discussed with emphasis on the modal 
decomposition of acoustic fields generated by PE models. 
Chapter III, Numerical Implementation, outlines the 


numerical modeling techniques used in implementing the c,- 


insensitive and double-precision WAPE models. 

Chapter IV, Numerical Results, explores the results 
obtained from each of the implementations examined. These 
include SPE, WAPE (single and double precision), and COIPE. 
Each model is examined for validity and accuracy against a 
number of standard test cases. The models were also analyzed 
for their ability to reproduce the standard normal mode 


basis set. 


das 


Chapter V, Conclusions and Summary, presents a summa- 
tion of conclusions and identifies areas of further investi- 


gation and research. 
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II. THEORETICAL BACKGROUND 


This chapter outlines the theoretical background sup- 


porting the development of the standard, wide-angle, and cj- 


insensitive implementations of the parabolic equation (PE) 
approximation to the acoustic wave equation. The bulk of 
this theoretical development is directed toward its applica- 
tion in the Monterey-Miami PE (MMPE) implementation of the 
wide-angle parabolic approximation. The MMPE implementation 
is emphasized here since it has been subjected to extensive 
analysis and numerous numerical improvements over its life 
span. It also serves as the framework of the COIPE model 
discussed and analyzed later. A summary of the method of 
normal modes with an emphasis on the modal decomposition of 
acoustic fields generated by PE models is also included. 
The entry point in deriving the PE is the Helmholtz 


equation, in cylindrical coordinates, 


2 2 
lop 2) (hop ep . pups 


Where D(r,z) gis the acoustic pressure, @kje= 7c, 1s the ref- 
erence wavenumber, n(r,Z,@) = C9/c(r,z,@) is the acoustic index 


of refraction, cg is the reference sound speed, and c(r, z, ) 


aE) 


is the acoustic sound speed. All of the environmental char- 
acteristics are represented in n(r,z,@). The treatment of 
density contrasts will not be presented here, but is 
included in the MMPE model by an effective index of refrac- 
tion term (Smith, 2000). A time-harmonic component of the 


acoustic field then has the form 


P(r, Zz, @, Mt) = p(r, z, p)e7 . (2023) 


A. OPERATOR NOTATION 


In order to develop the various parabolic approxima- 
tions to the Helmholtz equation, it is useful to introduce an 


operator notation. Let 


jp 2 (2.3) 
and 
OFn = (U+etvs ly, (2.4) 
1 a ie. 
where €=n?-1, = Gay, and v = 55.5. To account for 
kgdz kgr-0@ 


the dominate cylindrical spreading and to simplify the form 
of the Helmholtz equation, the acoustic pressure may be 


~a(r,z), which yields 
r 


alt 


defined as p(r,Z) = 


14 


(Pop + k§Qs,)a = 0. (a5) 
Proper factorization of this equation is accomplished by 


Getanang c= Gnu (Tappert, 1977), such that 


es + ikpQop)(Pop ~ ikpQ,p)u + iko[P,,, Qoplu =0Q0. (26) 
The commutator [ Pop: Qop] is exactly zero for layered media, 
and is assumed negligible. The remaining factors of Eq. 
(2.6) define the wave equations for the incoming and outgo- 
ing fields, the latter of which then satisfies 
Popul = IkpQop¥ (2e 7) 
or 


akg! = Oy ge (2.8) 


In cases where backscattered energy may be neglected, Eq. 
(2.8) provides the full description of the forward propagat- 


ing acoustic energy. 
B. THE SPLIT-STEP FOURIER ALGORITHM 


The most common methods of computing PE solutions are 
(1) the finite element method (Collins, 1994), (2) the 
implicit finite difference method (Lee, et. al., 1981), and 


(3) the split-step Fourier (PE/SSF) method (Hardin and Tap- 


ED 


pert, 1973). Of these three methods, the PE/SSF method has 
the distinct advantage of speed. This is accomplished by 
decomposing the acoustic field into a slowly modulating 
envelope function and a phase term which oscillates with the 
ikor 


acoustic frequency. Defining u = Ve and substituting into 


Ege 2.8) yields 


where H,, = 1-Q,, is a Hamiltonian-like operator which 


defines the evolution of the PE field function in range. 
In order to compute the solution to the acoustic field, 
a marching algorithm must be developed of the form 
Y(r+Ar) = O(r) V(r). (2.10) 
The propagator, @(r), is designed to propagate the solution 


out in range. The PE/SSF method accomplishes this via 
P(r) = e-ikpHop(r)Ar (2 : 1) 


where 


= 1 (r+ Ar) 
Hop(r) = | H,p(r)dr’ . (2p 
r 


In practice, the Hamiltonian is usually considered constant 


over a small range-step such that Hop(r) = Hop (r) - 
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The operator H,, (and Q,,) is not a scalar operator but 
a combination of scalar and differential operators. Each 
individual operator within H,, can, however, be applied by a 
Simple scalar multiplication in the proper domain. There- 
fore, each of the terms of H,, must be separated in order to 
comply with the PE/SSF algorithm. This necessitates the 


approximation of the square-root operator, specifically 


a 10° _. (a 1a 
Hop( so5 ane: n(r, 2)) = Tool 555) F U,, (ns, Z)) Tt Ceara ° 


(2913) 
It is this approximation that forms the basis of numerous 
higher-order PE model implementations. 


If the uncoupled azimuth approximation is employed 


then Vo, = 0 identically. The U,, is a multiplication opera- 


tor in z-space and is therefore a diagonal matrix. The oper- 


ator T,,. is not diagonal in z-space, resulting in different 
eigenfunctions being coupled. The Top is, however, diagonal 
in vertical wavenumber space. It is advantageous then to 
treat each operator separately, one in k,-space, and the 
other in z-space. Employing the Baker-Campbell-Hausdorff 
expansion (Bellman, 1964), 


eA+B = eAgBelA,B]+ A, [A, B]] + [B, [B, A]]} +... (2.14) 


wey 


where A = -ikjArT,, and B = -ikpArU,,. Top and Upp are both 
small, therefore it is assumed their products of second 


order are negligible. Using Eq. (2.14), ®(r) becomes 


4 OF oneal 
®(r) = cans 5) Uae + o0) ~ikoAuiecwaal 5 WAGs) 


(2g by 
Third order accuracy results from this centered step scheme 
(Jensen, et. al., 1994). Due to the formulation of the prop- 
agator there are no intrinsic losses due to the numerical 
scheme and energy conservation is retained (Tappert, 1995). 


In summary, the PE/SSF algorithm accepts the PE field 


function WY specified at some range r in the z-domain. This 


| — “ik Uo9(1) 
is followed by a multiplication of e 2 , the z-space 


operator defined at the beginning of the range-step. This 


product is then transformed to the k,-domain and multiplied 


by eo thoArT (k.) the k,-space operator. Another transformation 


to the z-domain and multiplication by the z-space operator 


~iky@ZU,,(1 + Ar) . 
e ~ , defined at the end of the range-step, produces 


the final resultant field function at r+Ar in the z-domain. 
Note that the FFT convention used here and in the numerical 


implementations is 


18 


c= FERCu(k) (2.16) 
Y(k,) = IFFT(‘Y(z)) 


Wrapping these steps into a single equation results in 


Wr+Ar,z) = (23 17) 


at . ike 
ED ea Ar, Z) . FET eto x ce Z) x V(r, |} : 


Cc. STANDARD PE DEVELOPMENT 


The lowest order approximation to the Hamiltonian oper- 
ator is commonly referred to as the standard parabolic equa- 
tion (SPE). This approximation is obtained by assuming the 
operators are small, 

€«l andyu«l : (27 18) 
The first of these is merely a recognition that sound speeds 
vary by less than 28 in most ocean environments, so n(r,z)=1 
everywhere. The second assumption can be seen to imply that 


kz 


vertical angles = snto are small. This limits the SPE 
0 


accuracy to small angles. 
With these approximations, the SPE is obtained by a 


binomial expansion of the square root operator, 


Ns 


1 l 1 1 
eee = 1-(1-5u-5e) = -Su-5e. (2. Dep 


The separated operators then have the form 


Q? 1 3° 
? | 3a 
wel ea) 2kgdz? a 
and 
Une z) = =3(n2r, Z)= 1) (23255) 


Note that in the context of Hj. acting as a Hamiltonian 
operator, the operators T,,. and Usp. correspond to kinetic 


and potential energy operators, respectively. In the verti- 


cal wavenumber domain, 


Tspe(k,) = 712 5 (2.220) 


D. WIDE-ANGLE PE DEVELOPMENT 


A higher order approximation to the operators was 
introduced by Thompson and Chapman (1983) and is often 
referred to as the wide-angle PE (WAPE) approximation. The 


operators of the WAPE are defined by 


Q° eee] 
Tyane(552) = 1-[+ eee Coen 
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and 
Sane: 7) —e(n(t, ai 1). (Ze24) 
In wavenumber space, the WAPE kinetic energy operator takes 


the form of 
‘ k,\241/2 
Twape(k,) = is[i-(= : ras) 


Note that energy propagating at angles beyond vertical, 


k,>kg, is evanescent since 


n k,\? 1/2 
Twape(k, > ky) = il) =1| (2.26) 


The wide-angle PE, and other higher-order approxima- 
tions to the SPE, was originally developed in an effort to 
extend the region of validity of the SPE. While the wide- 
angle PE succeeded in this goal, it introduced a number of 
other issues. These include a sensitivity to the selection 
of the reference sound speed and an inability to reproduce 
the standard normal mode basis set. This said, however, the 
wide-angle approximation has been very successful in produc- 
ing useful and valid results in many situations. The follow- 
ing section addresses an attempt to mitigate the reference 
sound speed sensitivity experienced by the WAPE approxima- 


EON . 


21 


E. COIPE DEVELOPMENT 


WAPE and other higher-order models are not exact repre- 
sentations and may under certain circumstances produce 
results worse than the SPE which they are meant to improve 
upon (Porter and Jensen, 1993). In these cases, it is often 
found that the error results from errors in the underlying 
normal mode basis set, their associated phase speeds, anda 


sensitivity to the choice of reference sound speed, cy. The 
choice of cp is the one ambiguous feature of all PE models. 


Due to this ambiguity, it is desirable to develop a model 
that is insensitive to the input reference sound speed. 


Gien,sany Significant deviations in the choice of cy would 


not adversely affect the computed results. In an attempt to 
develop a wide angle implementation that is insensitive to 
Go fapreriemerss Wyedeveloped the cpo-insensitive PE (COEPE) 
approximation. 

The COIPE approximation based on the WAPE approximation 


with the operator replacement 


a) 
ko(r)0z 


aca wig Z) 
kgdz? Ko(r)dz i 








(22 00) 
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where k)(r) is now range-dependent due to the range-dependent 
nature of the reference sound speed, c){r). Introducing the 


operator 





(2226) 


the cy-insensitive approximation to the Hamiltonian operator 


becomes 


Hys( Z) = (1-n(r,z)) + (1-1 -p(r)n7|(r, z)p(r)) . (2729) 


In this form the spherical wavefront property of the TC 
approximation is maintained for steep angles. Furthermore, 


Else approximation 1S not Sensitivesto the cholecsorac, and 


reduces to a c,o-independent approximation for small angles 


(Tappert, 1977). The resulting wave equation approximation 
is 

ov : 

a —ikgH,,,(r)'¥ . (27730) 


To implement this approach, a transform to a new depth 
coordinate 1s necessary. The transformation relationships 


are 


Zz 
Zz = | ac r)dr; (251) 
0 
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n(z,r) = n(z,¥n), (2320 


and 


W(z,r) = n-!/4(z, r)W(z, 18). (2338 


The corresponding Hamiltonian operator then becomes 


H(r) = (1—-n(z,r))+(1-1-p(r)) (2234) 
Note that in the tilde-domain this has the same form as the 
TC operator. The resulting form of the parabolic wave equa- 


tion in the tilde domain is 


aM = ik Hins(r)® - (2.35) 


The corresponding kinetic and potential energy operators of 


the COIPE in the tilde-domain may then be represented by 





T , K- =1- ;/l- ey 2.66 
Ins(1r . Kaa ( ) 

and 
Uins(r, Z) = 1—n(r,z) - (2 3a 


The COIPE requires, in general, the computation of co 


bersecach range step. To compute Cy, the condition that the 


water depth be unchanged in the z coordinate must be satis- 


fied. This is accomplished via 
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h(r) 
h(r) = | ain(zandz,, 25.0) 


0 
where h(r) is the water depth at range r. Under this condi- 
tion, the reference sound speed is 


h(r) —2 


a OS gg cies 
Co(r) = ao) SS” (Ze 39) 
0 


The University of Miami PE/SSF acoustic model (UMPE) 
has been tested using the technique outlined above with 
excellent results (Tappert, 1995). A numerical implementa- 
tion of the COIPE within the existing framework of the MMPE 
model is covered in the following chapter. Comparative 
analysis of this approach with the SPE, MMPE, and selected 
normal mode test cases is covered in Chapter IV. 

F. NORMAL MODE FUNCTIONS AND THEIR DECOMPOSITION FOR THE 
PE APPROXIMATION 

One of the difficulties encountered with the WAPE 
approximation is that it does not decompose into the stan- 
dard normal mode basis set as defined by normal mode theory 
and the SPE. Thomson (1993) looked at this problem for the 


SPE approximation. He showed that the field p that satisfies 
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the acoustic wave equation exactly relates to the field W 


satisfying the SPE approximation. The wavenumber of YW and 
the modal amplitudes may be determined exactly. The corre- 
sponding normal mode amplitudes and wavenumbers are then 
obtained. The WAPE, however, has modal eigenfunctions and 
eigenvalues of WY distinct from the acoustic wave equation. 

Smith and Smith (1997) showed that the modal decomposi- 
tion of solutions based on the WAPE generated erroneously 
range-dependent mode amplitudes in a range-independent envi- 
ronment. This was due to the mismatch between the standard 
normal mode basis set (based on the Helmholtz equation) and 
the proper basis set of the WAPE. Furthermore, they were 
able to develop an approximate numerical scheme for comput- 
ing the proper WAPE basis set. This resulted in the proper 
range-independent character of the mode amplitudes. Addi- 
tionally, they suggested that the higher-order COIPE version 
of the WAPE might overcome such issues and decompose prop- 
erly into the standard normal mode basis set. This will be 
one of the primary tests of the improved accuracy of the 
COIPE over the WAPE. 

Normal modes are obtained by using the method of sepa- 


ration of variables. As presented here, it assumes the ocean 
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is horizontally stratified, is of constant water density, is 
azimuthally symmetric, and is range-independent. As was 


presented for the PE, the time harmonic acoustic field 
p(t, zje7® at r>0O due to a point source at z = z, and r=0 sat- 


isfies the 2D Helmholtz equation (no azimuthal variations 


possible) 


2 
LO aD nO eee 
Of rsP) 4 SB + in (Z) p= 0 (2.40) 


For solutions of the form p(r,z) = 92) the depth- 
r 


dependent modal equation 
Z 
d 
Av (2) + (K§02(2) - K2)vq(2) = 0 (2.4) 


is obtained. Here, K2, 


the square of the horizontal wave- 
number for mode m, is the separation constant (eigenvalue) 
and V(Z) is®the specitic mode function (eigentumetion) 
associated with the separation constant. Given an arbitrary 


pressure field, the sum of the normal modes is required and 


is expressed as 


oo 


jalGe 74) + ial ; (2542) 


Za 


where the cylindrical spreading factor has been explicitly 

included. The factor 9,(r) can easily be shown to have the 

form of a Hankel function scaled by mode excitation values 

depending on the source strength and depth in the waveguide. 

In the far-field, the only range-dependence of this function 
is in the phase (having already accounted for cylindrical 


spreading). Thus, a generalized mode amplitude, A,(r), may 


be defined according to 


pu, z) = pu haleral® : (2.43) 


The normalized eigenfunctions form a complete, orthogo- 
nal basis set, defined by 


D 
Say Yal2 V2) a (2.44) 


0 
Thus, by using the orthogonality relation and assuming 
normalized modes, the generalized modal amplitudes can be 
extracted from the computed pressure field by forming an 
inner-product 


DAG) = afer = m(2) a (2.45) 
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This technique is used in the analysis and comparisons of the 
normal mode amplitudes of the SPE, WAPE, and COIPE models in 


Chapter IV. 


Finally, the corresponding depth separated equations of 
the various parabolic approximations introduced should be 


examined. First, note that the normal mode equation, Eq. 


(2.41), may be written in terms of the operators wu and € 


(defined in Eq. (2.4)) as 
Z 
(U+E)V_, = (S2- 1}: (2.46) 
0 


minis iundieates that the basis ser, V.;, are eigentuneurona oO: 


the combined operator (w+€) with eigenfunctions (K?,,/k29-1). 


The PE approximation has a depth separated equation of 


the form 
Kee 
opm = "Ky *™ r C2479) 


where the functions XY, are the eigenfunctions of the corre- 
sponding PE approximation. In order for the eigenfunctions 
1, and v,, to be the same (up to some normalization factors), 


the operator los must be some rational function of the oper- 


ator (u+é). For example, a simple polynomial series with 


Zo 


terms (w+e)", where n is some integer, will have the same 
basts=set Vv . 


The standard parabolic equation is a simple first order 


binomial expansion of the Hamiltonian operator, i.e. 
a 1 
oc ean +). (2 2a33 


This obviously satisfies the above criterion, therefore the 
SPE approximation has the same normal mode basis set as the 
fundamental Helmholtz equation. 

The wide-angle PE, on the other hand, may be written in 


the form 


Hyape = 2-V1+Hh- We (2.49) 
Performing a power series expansion on this operator to sec- 


ond order yields 
1 1 ia 
Hyape ~~ 5(U + €) + OCH + E*) (2:25:07) 


This shows that the WAPE approximation is only first order 
accurate in the operator (uw+€é) but has errors in the second 
order terms. It is ane weakness of the WAPE which generates 
some of the associated problems with this approximation. 


The cp-insensitive PE approximation, however, may be 


written as 
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He = 2 ae eee (2.51) 


Nile 


A power series expansion of this operator to third order 





yields 


el 1 2_ 1 3, We 
Hins = sb +E) + o(U +) g(t + ©) + 16° (252) 


Thus, the COIPE is accurate to second order in (W+€) with 
leading error in the third order terms. This was first shown 
by Tappert and Brown (1996). It is this feature of the COIPE 
which makes it an attractive alternative to the WAPE approx- 
imation. As will be shown, this also generates solutions 
which may be decomposed GS the standard normal mode basis 


set without significant error. 


Syl 
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III. NUMERICAL IMPLEMENTATION 


This chapter outlines the numerical implementation of 
the PE models used in the analysis of the cp-insensitive PE 
approximation (COIPE) as well as a double-precision version 
of the MMPE. The MMPE (Smith, 2000) is an existing well-doc- 
umented version of the WAPE approximation. In the cases 


where a SPE solution was required for comparison purposes, 


wn nn 


the operators Uyape and Twape were replaced with Usp. and Tspe 


within the framework of the MMPE model. As was discussed in 
the previous chapter, the COIPE implementation strives to 
eliminate the sensitivity to the chosen reference sound 
speed suffered by the WAPE approximation. As will be shown 
in the next chapter, a small change in the specified refer- 
ence sound speed can result in a dramatic change in the pre- 
dicted acoustic field. This error often increases with 
range and for complex, highly range-dependent environments 
can quickly result in the produced data being unusable. The 
numerical implementation described here was chosen such that 
it could be easily integrated into the existing MMPE model to 


provide the most reliable comparison of results. The 
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results of these comparisons are discussed in the following 


chapter. 


A. ENVIRONMENTAL GRID SIZES 


The MMPE, as well as the COIPE, implementation dis- 
cussed here uses a discretized grid to describe the environ- 
ment. The grid step-sizes are user definable and can 
dramatically affect the solution obtained by the model. The 
grid sizing method described here is that used by the MMPE 
model as well as the COIPE implementation of the MMPE. The 
environment is discretized by a mesh of size (Ar,Az). This 
results in the field and propagator functions being dis- 
cretized in depth with array length N, where N is the corre- 
sponding length of the FFT used in the PE/SSF algorithm. In 
other words, 

V(r, z) => P(r, Z,), (Sa) 
where 


res 1=1..i. , (3992p 


and 
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—— (3.3) 


The depth mesh is defined such that grid points lie on frac- 
tional values of Az. This convention results in the avoid- 


ance of carrying the zero-pressure value at z=0 throughout 
the calculations. It should be noted that because of using 
the full Fourier transform, half of the depth mesh values 
define an image ocean for negative depths. This has the 
added benefit of enforcing the surface boundary condition, 
defined in the next section, via FFT symmetry. Analysis by 
Tappert (1977) and Smith (2000) showed that by considering a 
lower limit on the allowable angles of propagation, a 
default value of Az may be defined and therefore a default 
transform size N. This is possible since the depth mesh 


influences the wavenumber increments Ak, via the FFT. Fora 


given environment with a relatively small angle of propaga- 
tion, the mesh size may be increased. An increase in the 
mesh size results in a corresponding decrease in the compu- 
tational time. As propagation angles increase, the mesh 


size is reduced. The most accurate solution, therefore, 
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should be obtained when both Az and Ar are on the order of a 


wavelength. 


B. TILDE TRANSFORMATION AND C, GENERATION 


The tilde transformation begins by defining a rescaled 


depth variable z in terms of the true depth z 


Zz 
z= [ae r)dz'. (3.4) 


0 
If it is required that the bottom depth remains fixed, in 


other words that 2 =z, =h, then the range-dependent bottom 


depth is 
h(r) h(r) h(r) 
Car 
h(r) = | A@ He = Soc ee' = eo | ee (3.5) 
0 0 0 
where Cy has been specifically declared as range-dependent, 
gue clchereromanaimor C(z;r). Solvangrior ¢p) (xr), 
h(r) —2 


_ | ae 
Co(r) = so tess (3.6) 
0 


This expression defines the rang-dependent reference sound 


speed to be used in the parabolic approximation. Its decla- 
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ration then eliminates the need for a user-defined value. 


Furthermore, since the definition is based on an evaluation 


of the local environment at each range, 


proc ty . 


If c(z) varies linearly with depth, 
c(Z) = c(z;_,)+b(z—z,;_,) = a+bz, 


where b is the sound speed gradient, 


Bee) ecem) 


— , 


Z—-Zi_] 
and 


then the integral may be written 


Z;- 


i-] 


Using a table of integrals, this becomes 


Z; 
Ja+bz’ 


=i 





Substituting into Eq. (3.11) for a and 5b, 





: 1 = i = ley a 
Ca . =Cj-1) (e263) (Nim i=) = 


7 


= =[ Jat bz; - Ja+bz,_1). 


it will preserve rec- 


(3.7) 


(3.8) 


eo) 


(3. £0) 


Cor eis) 


(oo) 


Now consider the sampling of the environment within the 


model. The typical definitions are 


2 = (k-5|Az p (3.733 
N 
k = 1,2,3...5 


where N is the FFT size. Based on this grid sampling of the 
environment, it would then be unusual if z,,, =h(r). Rather, 
the bottom lies just below the index 


nbi = int( X22) F (3.14) 
Az 


such that z,,sh(r). To perform the integrals needed, how- 


ever, the conditions 


2 = We (32259) 
ZNs = h(r) , 


must be satisfied. This may be accomplished by first trans- 


forming to a Zz; space (which is not equally spaced), which 


ealistres Equa. 15), andethenminterpolating zjeto z, where 
3 l 
2 = (k-5)d2. (35g, 


Part of the difficulty is to define the sound speed on 
the non-grid points z=0 and z=h(r). Fortunately, this is 


essentially already accomplished in the MMPE model. The 
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surface sound speed is required as an input to the MMPE, and 
is stored in the sound speed array by the environmental prop- 
agator subroutine. The sound speed at the bottom depth, 
h(r), is also already estimated at the grid point nbi. While 
this is an estimate, it provides a good approximation for 
cases where energy doesn’t significantly interact with the 


boctom. This results in 


Zz} = Oz = 0=>c(0) = SS , (3257) 
Zn, = h(r)>z = h(r)>c(h(r)) = SSBW, 


where SSO is the surface sound speed input into MMPE and SSBW 
is the estimated water sound speed at the water/bottom 
interface. 

The numerical scheme for computing cy and the tilde 
transformation can now be defined. Using Eqs. (3.6) and 


wo 2), 


Co(r) = (3.96) 


: -2 
nbi 


call = |: y Buel -| h(t) ~ Znpi | 
Je, +VSSO} | Jep+ Jey 1} LVSSBW + Jenpi 
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where it is assumed that the sound speed has already been 


interpolated onto the gridded mesh z, at values c,. Simi- 


larly, Eq. (3.4) with z,; = 0 becomes 


z2 = 2,/co, (3.19) 
— Z.—Z 
zj = 2)-1 + 2,feg——“E4, j= 3...nbit+1, 


ZNs = Zns-1 + 2,/¢g ——>—— 
JSSBW + Hes 


Note that Ns = nbi+2 since the points z=0 and z=h have 
been added. As a result of this transformation and the above 


definitions, 


c, = SS, (3.20) 
CNs = SSBW, 


Cj 1, NsiemC ey 
Given these relationships, it is now necessary to 
interpolate in depth the tilde domain sound speed onto the 


standard grid. Defining 


Zz, = (k-5)az. (3.29 
N 
Kes. 2: 35 
it is desired that 
C(Zk) = C(Z,). (3222) 


40 


A simple linear interpolation from c(zj) to c(z,) could be per- 


formed. This, however, would not preserve the linear gradi- 
PHe ndeiise Orec(2)), end themrelatrronedeeimed ss yng. (sma2) 
would not be satisfied. Instead, assume 

Zj-1 SZ <2, (328) 


and then define the fractional depth as 


ka Zi 
f=. (324) 
Zj— 2Zj-1 


Tappert (1995) shows that the inversion of the mapping 
needed produces 
C= Cai 2. (Cee (3.25) 


where 


fi, 
&, = Fela Se (3.26) 


Cc. ACOUSTIC FIELD GENERATION 


In the MMPE model the source amplitude is defined rela- 


tive to a small finite distance from the source. This is 


done to avoid the singularity at range r=0 in 


pS Po,|Sowe' (3.27) 
rT 
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Therefore, the source is defined such that p(r=Rp) = Pp. The 
reference range is defined to be one meter to remain consis- 


tent with most sonar equations. Using this convention, the 


source level is related to Py by 


P 
SL = 20l08( 5) dB re P. Ry (3.28) 


T 


where wr =——t Nba saet tne reLerence range OL Ry. 


The form of the source field V(r=0,z) must now be 


determined. Rewriting Eq. (3.27) as 


ikor — 1 [x 
Wr, z)e P, pee ee (35295) 
then 
: ae 
Y(r=0,z) = lim = /—p(r,z), 3.50 
(= 0,2) = lim p- fe-pC. 2) (3.30) 


where in the vicinity of the source the pressure filed takes 


the form of the spherical Green’s function. If |p| = re ne 
required at R=Rpj), then let a= Pp9Ry and therefore 
— —4_pikyR 2am 
P4nR Se 


where R = jr*+z*. The source is represented as a point 


source at (0,z.) by 
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V(r =0,z) = ad(z-z,). (3. 


Integrating Eq. (3.30) over all depths of the real ocean 


results in 


= ae eikor 
ms NR (Re ane 7a ale i 


In the far-field, r»z, R may be approximated such that 


re 


Phus reducing Eq. (3.33) to 


oO 





a= /Ro lim ~ ee “Fide ~ {Ro lim Seili , (3. 
r>04n./r lim wiles 


0 


and leaving 


_ Rg 
rod : (3 
27ky 


Se ) 


rosy) 


wo 4) 


58)) 


236) 


Performing a Fourier transform of Eq. (3.32) with the 


influence of the image source included, and specifying the 


source in the k,-domain yields 


‘ir = 0k) = af [5(z-z,) -8(z + z,)Je7tk:"dz = —2iasin(k,z,) . (3 


OO 
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ome 


Equation (3.37) states that the starting field, in the wave- 
number representation, has constant amplitude moduua Med by a 
phase due to the interaction of the source and its image. 
This representation is consistent with an omnidirectional 
source that supplies equal energy to all wavenumbers. 

Rather than setting the amplitude of this function to 
unity to ensure equal population of all directions, a smooth 
taper is used at high absolute wavenumber values to limit the 
angular dimension of the source and to reduce sidelobe 


influences. This is necessary due to the wide-angle approx- 
imation being valid to only approximately 40°, and due to the 
Pesesecci1on on how large the finite FFT will allow K,/Kg- 
The wide-angle source is modified to produce an accurate 
solution in the far-field (Thomson and Bohun, 1988) by 
2\-1/4 

Eke) = 1-3 5S Slay (3.38) 

The limits on k, are chosen as such since for |k,|>kg the 


angles of propagation are imaginary. Therefore, the source 


kK, <li 


function is tapered within the limits . 
0 
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Finally, to account for the half-mesh symmetry of the 


depth grid a phase term in the wavenumber domain of e “2 is 
added. The wavenumber domain starting field for the wide- 


angle point source can then be written as 





‘ iR 2\- 14 5 Az 
W(r = 0,k,) = -2i Fat sink 21-73) ee: (3.39) 


Following Tappert’s procedure, for the COIPE it is eas- 
iest to generate the starting field in the tilde domain. The 


starting depth, z,, must then be transformed to tilde space. 


Defining 


aia) (3.40) 


the source depth in the tilde-domain is then 


~ ~ Zi) — Zi_ 
Zs = Zj-1+ (fe; + fo;_,) f= 3041) 


Ju; + Je | 
where 
Note that Eq. (3.42) performs a linear interpolation of c at 


the source depth. Since, from Eq. (3.22) 
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then Eq. (3.41) becomes 


~ ~ 


Zp = Zj-1 


or 


which is merely Eq. 


2/e(a-7-1)f _ - 


4 2aleol2s- 25-1) 


fe; + fej—1 


(3.19). rewr Ebenmuicingec. = C(Z.).. 


(3.433 


(3.44) 


(3.45) 


Also 


from the expression for the tilde domain transformation, 


it ELollows that 


Therefore, 


So 


Z = face 


0 
= =a Zan) 


dz 
=~ = (n(z,r))-!/4 
dz 


= Jaw, r))-1/2¢7’ . 


0 
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(3.46) 


(3.47) 


(3.48) 


(3.49) 


Since c(z) = c(z), then n(z) = n(z), such that 


~ ~ 


Zz Zz 
z= fone. ae = [ deceyaz (3.50) 
0 





NCo(T) 


0 
Equation (3.50) could be solved numerically, but it does not 
have the analytical expression the previous transformation 
did. Rather, the previous transformation for source depth 


can be applied to receiver depths at the original gridding of 


the environment. In other words, consider 
1 
Z, = (k-3)az C3ie5 1) 


which is essentially what was obtained in Eq. (3.19), where 


now the tilde depths are those found in the original trans- 


Formation vector, z; for j=2...Ns—l1. Below Zys-1, Ze2Nns-1 15 


used. 


To obtain the field, however, it must be noted that 
W(z,r) = [n(z, r)]-!/4¥ (z, r) (3.52) 
Or, Since n(z) = n(z), 
[n(z, 14 w(z, r) = V(r). (3.53) 
In either PPR (3) must first be interpolated from W(z,) - 


Since there is no simple gradient structure for Y(z), a 
ay 


numerical scheme must be used. For this case, a typical 5- 
point interpolation scheme is used to obtain (z;) . Li Eon 


(3.52) is used, then n(z,,r) is evaluated and the values of 


(z,,7) = [n(Z,, 1) Y(z;) (3.54) 


are computed where k=j-1 and c(z,) is evaluated from the 


original profile. Note that [n(z,,r)]4¥(z) could have been 
computed and then interpolated to 2}. However, this would 


only be an approximate interpolation of n. Therefore, the 


former method is an exact interpolation of n and only 


approximate in W, which is expected to be the more accurate 
method. 

As was briefly discussed in the previous chapter, the 
COIPE approach is expected to remove the reference sound 
speed ambiguity present in the WAPE. Removing this sensi- 
tivity should provide more accurate results without the need 
to determine an appropriate reference sound speed for the 
specific problem. In the following chapter the COIPE is 


benchmarked against the SPE, WAPE, and normal mode theory. 
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IV. NUMERICAL RESULTS 


The previous chapters present an outline of the histor- 
ical and theoretical background of the parabolic equation 
approximation to the acoustic wave equation as well as the 
COIPE implementation. In this chapter, the focus shifts to 
the analysis of attempts to improve the accuracy and valid- 
ity of the WAPE approximation. First, the effects of improv- 
ing the numerical accuracy of the WAPE, as implemented by the 
MMPE model, is evaluated. Second, the COIPE implementation 
is evaluated under a number of test cases. Finally, the SPE, 
MMPE, and COIPE are subjected to the normal mode decomposi- 


tion technique previously discussed. 


A. DOUBLE-PRECISION IMPLEMENTATION OF THE MMPE 


While the MMPE has been used successfully under a wide 
variety of test cases and environments for a number of years, 
it was recently tested extensively on a number of specific 
problems. These tests were conducted as part of the Shallow 
Water Acoustic Modeling (SWAM) Workshop held at the Naval 
Postgraduate School in 1999. The SWAM workshop provided a 


set of detailed environmental cases that served as the basis 
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for comparison between the numerous modeling efforts being 
undertaken in the United States and internationally. As 
part of this workshop, the convergence sensitivity of the 
MMPE to the choice of the range step, Ar, was analyzed 
(Smith, 2000). As a result of this analysis, which is dis- 
cussed in more detail later in this section, an effort to 
develop a more numerically precise MMPE implementation was 
undertaken. Increasing the numerical precision of the model 
would also help answer whether truncation of the field, as 
postulated in the introduction, is a source of error in the 
model . 

In an effort to avoid introducing new errors or anoma- 
lies into the MMPE, no algorithm changes were made to the 
program code itself. Rather, the variables and array decla- 
rations, as well as the corresponding post processing rou- 
tines, were recast as double-precision. The original MMPE 
code was, with the exception of the FFT subroutine, predomi- 
nantly single-precision. One of the concerns raised by this 
procedure was that computational times would increase dra- 
matically. This, however, was not the case. Since the bulk 
of the computational work is done by the FFT subroutine, 


which was already in double-precision format, computational 
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times only increased by three to five percent. These 
increases are deemed relatively insignificant and do not in 
and of themselves hinder the use of the double-precision 
MMPE (DP-MMPE). 

As a method of comparison with the SPE and existing 
MMPE, the DP-MMPE was tested using one of the SWAM test 


cases. The test case chosen was referred to as SWAM Flat-a. 


The environment has an isospeed water column, Cy = 1500 m/s 


Po = 1.00 g/cm?, overlying a flat bottom. The bottom proper- 


ties are defined every 2 km and are linearly interpolated out 
to a range of 20 km. Compressional attenuation is fixed at 
0.1dB/A and there is no shear in the bottom. The source is 
at a depth of 30 m and a frequency of 250 Hz CW. The envi- 
ronment is shown below in Fig. 4.1. Fig. 4.2 shows the full, 
CW field at 250 Hz for the source at 30 m as obtained from 
using the MMPE model solution. 

During the theoretical development of the wide-angle PE 
approximation it was noted that the range step size Ar, as 
based on the centered-step scheme, is roughly third order 
accurate. The implications of this fact are that as the 


range step size is decreased, the solution accuracy should 


Sal 


improve until a stable solution converges. 


models do follow this behavior. 
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Figure 4.1. 
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The Flat-a test case environment. 


To illustrate the relation between the range-step size 


and the covergence of the solution, 


were generated for the Flat-a case discussed above. 


a number of solutions 


The only 


parameter varied in any of the cases was the size of the 


range-step. 
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Figure 4.2. MMPE field solution for Flat-a case. 


for a receiver at 35 mis shown in Fig. 4.3. The upper and 
lower panels of Figs. 4.3 and 4.4 show the first and last 5 
km of the solution, respectively. Note that the solution 

appears to be converging as Ar~5 m. As can be seen in Fig. 
4.4, however, the solution diverges for further decrease of 


Ar; the solution degrades for values of Ar less than about 


5m. The solution appears to reach convergence at Ar~d. 
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Figure 4.3. Convergence testing for various range-steps. 
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Figure 4.4. Stability testing for various range-steps. 
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A possible cause of this divergence was postulated to 
be truncation of the field and a lack of numerical precision 
in the declaration of variables in the MMPE computer code. 
In an effort to test this theory the DP-MMPE model was run 
with identical model inputs to those used in the convergence 
and stability testing shown above for the MMPE. The results 


of this testing are shown below in Figs. 4.5 and 4.6. 
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Figure 4.5. DP-MMPE convergence testing for various range- 
step sizes. 
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Figure 4.6. DP-MMPE stability testing for various range-step 
sizes. 


Clearly, there does not appear to be any discernible 
improvement in the convergence of the solution through the 
use of the DP-MMPE. It exhibits the same general behavior as 
that of the original MMPE model. A comparison between the 
MMPE and DP-MMPE at the 5 m and 2 m range-step sizes is shown 


meow in Fig. 4.7. 
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Figure 4.7. DP-MMPE and MMPE comparison for range-step sizes 
of 5 mM and 2 m- 


Referring to Fig. 4.7, it can be seen that for a range 
step size of 5 m the solution produced by the two models are 
indeed overlapping. For a range-step size of 2 m the two 
solutions differ slightly with range, with the DP-MMPE 
appearing to diverge less at some range points. One may con- 


clude, however, that while the DP-MMPE does no worse, it does 
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not offer a solution as to why the PE/SSF algorithm seems to 
have a lower limit on the size of the range-step. 

As a result of separate analysis (Smith, 2000), the 
primary cause of this non-convergence was found to be 
related to the structure of the propagator functions. 
Recall that the environmental propagator function is of the 


-ikpArU,,(r, z) 


£Orm € and the wavenumber propagator function is of 


~ikpArT,,(k,) 


the forme The wavenumber propagator decays expo- 


nentially beyond k,>ky. From this analysis it can be shown 
that the PE/SSF algorithm, for large range-step size, 
attempts to include a large amount of phase information into 
each range-step. If Ar is chosen to be too large, then 


errors are generated in the solution. Conversely, as Ar is 
reduced in size there is inadequate phase information con- 
tained in each range-step and the solution once again 
becomes inaccurate. The most accurate solutions may be 
expected for range-steps where a full cycle of phase infor- 
mation is contained in each propagator. The level at which 


Ar appears to provide the greatest degree of convergence is 


Ar~X (Smith 2000) for shallow water environments such as 


tested here. 


Sis: 





Before leaving the area of solution convergence based 
on the range-step size, a final comparison between the MMPE 
and the COIPE is made. While the COIPE model was not imple- 
mented specifically to address this problem, it is of inter- 
est as to whether any of the changes made in the COIPE affect 
the convergence due to Ar. Below, in Fig. 4.8, the COIPE 
model is compared against the MMPE for the same test case 
used above. The MMPE results were chosen over the DP-MMPE 
Since the DP-MMPE solutions did not offer any definitive 
advantages over the normal implementation. As Fig. 4.8 
shows, the COIPE does not appear to offer any greater degree 
of convergence for the given Ar sizes in this environment. 
These results are inconclusive in regards to the COIPE 
model, but the implementation was not expected to affect the 
Ar sensitivity. 

Finally, as a means of comparison against another mod- 
eling technique for this test case, the MMPE results for 
range-step size 5 m are compared to those developed by Mikhin 
(2000) for SWAM. These results were chosen as a benchmark 
due to the high degree of agreement with several test cases 


and several modeling techniques. 
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Figure 4.8. COIPE and MMPE convergence testing. 


The comparison of the MMPE solution with the solution 
provided by Mikhin is shown below in Fig. 4.9. Solutions for 
rang-step sizes of 20, 5, and 2 m were chosen to demonstrate 
the degree of convergence with the reference solution for 
varying Ar. The figure shows very good agreement between 


the MMPE and Mikhin solutions at short ranges. The two tech- 
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niques show a greater degree of separation with increasing 
range. This is most likely the result of cumulative phase 
errors inherent to the wide-angle PE technique employed 


here. The figure also demonstrates, once again, that the 
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Figure 4.9. MMPE and Mikhin solution comparison for various 
range step-sizes. 


choice of a 5 m range-step size provide the greatest degree 


of agreement with the reference solution. 
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B. THE COIPE IMPLEMENTATION OF THE MMPE 


As was discussed in the last chapter, the COIPE model 
was developed by Tappert (1995) in an attempt to remove or 
reduce the dependence of the wide-angle PE approximation on 
the choice of the reference sound speed value while increas- 
ing the overall accuracy of the model. While this sensitiv- 
ity has been observed in a number of test cases, it was 
extensively studied during one of the PE Workshop II (PE-ITI) 
test cases (Porter and Jensen, 1993). This case, referred to 
as the Porter duct problem, was based on the observation that 
many of the wide-angle split-step PE approximations overes- 
timated the value of the TL when applied to long-range prop- 
agation in a leaky surface duct environment. The SPE 
implementation of the PE/SSF algorithm does not generate 
this anomaly. Additionally, other PE implementations based 
on Pade’ series expansions of the square root operator 
(which maintain the proper normal mode basis set) do not 
result in this overestimation of the downrange TL. The 
results of PE-II showed that many wide-angle models that do 
not use the Thomson-Chapman (TC) wide-angle approximation 


avoid this problem. 
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The Porter duct environment is one of a source and 
receiver located in a 250 m deep surface duct. The source is 
at a depth of 25 m and frequency of 80 Hz. The receiver is 
at a depth of 100 m. The ocean bottom is range-independent 
with a lossy fluid half-space beginning at a depth of 4 km. 
Of the 78 propagating modes in the water column, one is 
trapped in the surface duct. Figures 4.10 and 4.11 illus- 
trate the sound speed profile and TL field for this environ- 


ment respectively. 
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Figure 4.10. Porter duct sound speed profile plot. 
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Frgm@resd4ei1. Ti, field plot for the Porter duct problem: 
Finn Jensen, during PE Workshop II, showed that the SPE 
was in excellent agreement with the normal mode reference 
solution used in the test case (Porter and Jensen, 1993). 
Because this reference solution is no longer available, the 
SPE results were used as the reference for this analysis. 
The only alterations to the various model input parameters 
changed in the following analysis is that of the reference 


sound speed, Cg. Prior to looking at the COIPE model or the 
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effects of varying cy, the TL plot for the receiver depth of 
100 m is shown in Fig. 4.12 for the SPE and WAPE implementa- 
tions when a reference sound speed of cy = 1500 m/s is used. 
It can be clearly seen that the WAPE solution shows a dis- 
Lt is@ehis 


tinct drop off at ranges beyond about 60 km. 


anomaly that the COIPE model is meant to eliminate. 
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TL Plots for SPE and WAPE at the receiver depth 
o£ LOC m andcy = 1500 m/s. 


PrGgure 4.12. 


The reference sound speed was then set at 1485 m/s and 


the two models were run again. These results are shown in 
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Pige4- b3eeethis plotgmillustrates the cttect™the valuctomec, 
has on the WAPE solution. By “correctly” choosing the value 
of Cp, the WAPE correctly determines the long range TL. The 
difficulty, however, is in correctly determining this value. 
A change in co by as little as 5 m/s 

to either side of this “correct” value results in the long 


range drop off shown in Fig. 4.12. 
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Figure 4.13. TL Plots for SPE and WAPE at the receiver depth 
of 100 m and cp = 1485 m/s. 
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The COIPE was then tested using identical environmental 
inputs to those used in the analysis above. COIPE model runs 
were made at Cp values of 1500 and 1485 to provide comparison 
with the SPE results. These are shown in Figs. 4.14 and 4.15 
respectively. The COIPE model compares well with the SPE 
solution at all ranges and there is no drop off at extended 
range aS was observed with the WAPE. The COIPE solution is 


indeed insensitive to changes in the reference sound 


TL (dB re 1m) 





Figure 4.14. TL Plots for SPE and COIPE at the receiver depth 
©fsl00°meand <c, = 1500) m/s. 
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Figure 4.15. TL Plots for SPE and COIPE at the receiver 
depth of 100 m and cy = 1485 m/s. 


speed, as shown in Fig. 4.15. While both plots in Fig. 4.15 
do show phase differences from those of Fig. 4.14, the COIPE 
follows the reference solution equally well in both cases. 
The COIPE was tested at a wide variety of cg values with sim- 
ilarly successful results. 

Clearly the COIPE corrects the TL drop-off encountered 


in the MMPE and similar wide-angle implementations. The 
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COIPE corrects the small error in the phase calculation 
caused by the incorrect choice of the reference sound speed. 
A small error in this phase calculation can cause a large, on 
the order of 20 dB, downrange TL error (Porter and Jensen, 
1993). The test case highlights this problem because it sup- 
ports a single trapped mode in the surface duct that is a 
leaky mode. The trapped mode continually loses energy into 
the lower region where, due to the strong upward refracting 
profile, it is directed back up toward the surface duct caus- 
ing interferences. Small phase errors result in incorrectly 
applying these interferences to the surface duct mode. The 
transformation and subsequent calculation of the reference 
sound speed at each range step removes the source of these 
phase errors and thus correctly applies the interferences 


described above. 


Cc. NORMAL MODE DECOMPOSITION OF THE PE FIELD 


The mode functions of the wide-angle PE approximation 
form a different basis set for modal expansion than those 
obtained using standard normal mode theory. As discussed in 
the theoretical development presented earlier, the SPE does 


reproduce the standard basis set. This section presents an 
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analysis and comparison of the modal decomposition of the 
SPE, WAPE, and COIPE models. 

The ocean environment chosen for this analysis is the 
Munk canonical profile (Munk, 1974). This is a deep ocean, 
range-independent environment with sound speed profile char- 


acterized by 
[c(z)-—c¢g]/cg = E(e1-H-1), (40) 


where € = 0.0057, the scaled depth variable n = 2(z-z 7B, 


ae) 
and the reference sound speed is 1490 m/s. The axis depth 


was set at Z = 1000m, and the bottom depth at 4500 m. The 


axis 
sound speed profile is shown below in Fig. 4.16. The source 
was placed at a depth of 1000 m with a frequency of 100 Hz. 
A computational range of 100 km was used for PE field calcu- 
lations. A bottom depth of 4500 m and a computational depth 
of 8000 m were used. 

The acoustic pressure field was calculated using each 
of the three models. The field is shown in Fig. 4.17. The 
resulting fields were then decomposed into normal modes and 


the mode amplitudes determined using the techniques outlined 


in Chapter II for each of the models. Since the 
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Figure 4.16. Munk canonical sound speed profile. 
environment chosen is range-independent, the mode amplitudes 
should remain constant with range. The amplitudes for 
selected modes are shown in Figs. 4.18 thru 4.20 for the SPE, 
WAPE, and COIPE respectively. Modes were chosen to avoid 
overlapping plots to improve readability and to avoid bottom 
interactions. The reference sound speed is 1500 m/s for this 
series of plots. Each of the three models produce satisfac- 


tory results for the lower modes. Note the nearly constant 
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Figure 4.17. TL Field plot for Munk canonical profile. 
amplitude with range for all three models. The higher order 
modes of the WAPE plot, however, show increasing levels of 
fluctuations. Conversely, the COIPE modal amplitudes com- 
pare well with the SPE results and do not show the high 
degree of fluctuation seen in the WAPE plots. 

The three models were then tested for this environment 
using a reference sound speed of 1490 m/s. The results are 
shown in Figs. 4.21 through 4.23. The results for the COIPE, 
Fig. 4.23, are visibly unchanged, as expected. The WAPE 
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Figure 4.18. SPE mode amplitudes with range for modes 1, 5, 
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Figure 4.19. WAPE mode amplitudes with range for modes 1, 5, 
Opals 20pm, SOgmson 40, 25> and 50. 9(c, = 1500 m/s) 
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Figure 4.20. COIPE mode amplitudes with range for modes 1, 5, 
UOjeeeore eo OU, 35, 40, 45, and 50. (cp = 1500 m/s) 
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Figure 4.21. SPE mode amplitudes with range for modes 1, 5, 
DOVerS, 2emeeco, 30, 35, 40,845, and 5028 (cy = 1490) m/s) 
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Figure 4.22. WAPE mode amplitudes with range for modes 1, 5, 
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Figure 4.23. COIPE mode amplitudes with range for modes l, 5, 
yee Opec os O poo oO. 45, and 50. (cy = 1490 m/s) 
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amplitudes show somewhat less fluctuation for the higher 
modes. This can most likely be attributed to the sensitivity 
to the reference sound speed discussed in the previous sec- 
ELone 

This analysis has demonstrated the wide-angle PE’s lack 
of the ability to reproduce the standard normal mode basis 
set as defined by normal mode theory and as duplicated by the 
SPE. The WAPE shows increasing levels of fluctuation in the 
modal amplitudes with increasing mode numbers. This fluctu- 
ation results in the produced PE field not exactly repre- 


senting the environment under study. The COIPE, however, 


does not show the same degree of fluctuation in modal ampli 
tudes, and therefore should better represent the acoustic 

field. It should be noted that the modal amplitudes of the 
COIPE do not exactly match those of the SPE. This is due to 
the generation of the source function in the tilde-domain. A 
more proper treatment might generate the starting field in 

real space and then transform to tilde space. However, this 
would introduce some numerical errors to the transformation 
of the filed form tilde to real space. This is an area for 


further analysis and study. 
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V. CONCLUSIONS AND SUMMARY 


This thesis has examined the parabolic equation approx- 
imation to the acoustic wave equation and, in particular, 
two possible implementations aimed at removing inherent WAPE 
errors were analyzed. The theoretical background underlying 
each of these models was presented along with the specific 
numerical implementation used for testing and comparison 
with existing PE approximations. Finally, numerical results 
were provided comparing the SPE, WAPE, DP-WAPE, and COIPE 
models against a number of test cases. These included 
results from the SWAM workshop, the Porter duct problem, and 
the Munk canonical profile. 

The parabolic equation and many of its derivatives suf- 
fer from a number of anomalies or errors depending on the 
particular implementation being used. Of these, three were 
examined in this analysis. First, the selection of the 
range-step size used by the PE/SSF algorithm affects both 
the convergence and the stability (Smith, 2000) of the solu- 
tion with an error to second or third order in the range-step 
increment (Jensen, et. al., 1994). Second, the WAPE suffers 


an anomaly that under certain environmental conditions it 
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incorrectly computes the down-range TL (Porter and Jensen, 
1993). This anomaly is manifest in a distinct drop in This 
Results of the Parabolic Equation Workshop II (PE-ITI) 
attribute this error to incorrectly handling the loss of 
energy from a leaky surface duct mode and its subsequent 
interaction with the remaining sound field. It was shown 
here, and in PE-II, that the manifestation of this anomaly is 
highly dependent on the choice of the reference sound speed. 


A small change in the value of cp) results in a significant 


change in the computed down-range TL. Third, the WAPE does 
not reproduce the standard normal mode basis set as defined 
by normal mode theory and as reproduced by the SPE. In addi- 
tion to the phase errors experienced to some degree by all PE 
approximations, the WAPE mode amplitudes exhibit a range 
dependence in a range-independent environment (Smith and 
Smith, 1997). This effect was found to contribute to the 
anomalous TL findings just described. 

A double-precision version of the WAPE model was imple- 
mented in an attempt to remove or improve the convergence 
sensitivity to the selection of the range-step size used by 
the PE/SSF algorithm. The double-precision model was tested 


against the existing SPE, WAPE, and COIPE models in a simple 
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range-dependent environment originally defined for SWAM. 
The double-precision model was found to do as well as, but no 
better than, the existing approximations under the tested 
conditions. From this analysis it was concluded that the 
range-step convergence sensitivity is inherent to the PE/SSF 
algorithm used by the models. This is supported by separate 
analysis done by Smith (2000). The most accurate solutions 
may be expected for range-steps where a full cycle of phase 


information is contained in each of the propagators used by 
the algorithm. The level at which Ar appears to provide the 


greatest degree of convergence is Ar~A (Smith 2000) for 
shallow water environments such as those tested. The COIPE 
model was also tested for sensitivity to the choice of the 
range-step on solution convergence. It was found to be sim- 
jlarly affected, as was expected. 

The COIPE model, a WAPE implementation, was developed 
to remove or reduce the solution sensitivity to the choice of 
the reference sound speed value (Tappert, 1991). The COIPE 
model uses an operator transformation which results ina 
range-dependent series of values for the reference sound 
speed. The COIPE model was compared with the SPE and WAPE 


models using the Porter duct problem developed for PE-II 
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(Porter and Jensen, 1993). These tests were conducted using 
a variety of reference sound speed values, with successful 
results in all cases. The COIPE model does not exhibit the 
TL drop with range seen with the WAPE. 

Finally, the COIPE was tested against the SPE and WAPE 
models for its ability to decompose properly into the stan- 
dard mode basis set. The models were tested in a deep-ocean 
Munk canonical profile environment. As discussed, the WAPE 
does not correctly decompose with this basis set, and shows a 
range dependence in its modal amplitudes in a range-indepen- 
dent environment. This range dependence can be seen as fluc- 
tuations in the modal amplitudes with range. The COIPE modal 
amplitudes show much less fluctuation with range than those 
of the WAPE, and more closely reproduce those obtained by the 
SPE. This is a result of the improved handling of the phase 
by the COIPE as outlined above and demonstrated through the 
modal decomposition. 

The double-precision model testing provides further 
evidence, in addition to that shown by Smith (2000), that the 
small range-step sensitivity is a fundamental property of 
the PE/SSF algorithm, and no further investigation of this 


property is necessary. The results of the COIPE model test- 
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ing, while promising, are not conclusive. The COIPE cor- 
rects the error seen in the Porter duct problem, and improves 
the reproduction of the standard normal mode basis set, but 
requires further investigation. Preliminary testing in 
shallow water environments indicate the COIPE may not pro- 
duce results significantly different from those of the WAPE 
model. Additionally, the mode amplitudes themselves, while 
showing less fluctuations, are of different relative magni- 
tudes from those produced by the SPE. This remains as an 


area for further investigation. 
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